Melting of Rare-Gas Crystals: Monte Carlo Simulation versus Experiments 
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We study the melting transition in crystals of rare gas Ar, Xe, and Kr by the use of extensive 
Monte Carlo simulations with the Lennard- Jones potential. The parameters of this potential have 
been deduced by Bernardes in 1958 from experiments of rare gas in the gaseous phase. It is amazing 
that the parameters of such a popular potential were not fully tested so far. In order to carry out 
precise tests, we have written a high-performance Monte Carlo program which allows in particular 
to take into account correctly the periodic boundary conditions to reduce surface effects and to 
reduce CPU time. Using the Bernardes parameters, we find that the melting temperature of several 
rare gas is from 13 to 20% higher than that obtained from experiments. We have throughout 
studied the case of Ar by examining both finite-size and cutoff-distance effects. In order to get 
a good agreement with the experimental melting temperature, we propose a modification of these 
parameters to describe better the melting of rare-gas crystals. 

PACS numbers: 34.20.Cf, 64.70. Dv, 65.20.+W, 65.40.-b: 



I. INTRODUCTION 

Melting of crystals has always been a fascinating sub- 
ject for more than a century since the discovery of the 
empirical Lindemann's criterion [l|. The Lindemann's 
criterion says that if the average of vibration amplitude 
u, namely y/ (u 2 ) , exceeds a certain value, usually 10% of 
the distance between nearest-neighbors, then the melting 
occurs. Times and over again, many authors have tried 
to find out microscopic precursor mechanisms that lead 
a crystal to melt. Until 30 years ago, one of the favorite 
pictures of melting is the softening of a phonon mode 
due to the temperature (T ). The atoms have no longer 
restoring forces which keep them staying close to their 
equilibrium positions: they move around and the system 
goes to a liquid state. The soft-mode picture has encoun- 
tered some scepticism because in real crystals as well as 
in simulations one observes that well below the melting 
temperature (T m ), many defects, dislocations, intersti- 
tial atoms ... are excited. Therefore, it is hard to believe 
that the system stays in a periodic structure with prop- 
agating phonon modes up to T m . Evidence of defects 
is found in many works [2]-|5(. Another question that is 
unsolved in a clear manner is the form of the potential 
that binds the atoms together in a given lattice struc- 
ture. In a microscopic point of view, the potential should 
come mainly from the symmetry of atomic orbitals. But 
ab-initio calculations arc still far away from being able 
to use realistic hypotheses @ . Empirical potentials have 
been used instead to study melting. One can mention the 
popular 6-12 power Lennard- Jones (LJ) potential @, 
various similar power potentials, the many-body Gupta's 
potential [§], the Stillinger- Weber (SW) potential [To(|. 



and the Tersoff potential |ll|, [l2[ . Two-body potentials 
such as the LJ one crystallize atoms in the FCC at low 
temperatures and nothing else; this comes from the fact 
that LJ potential is isotropic so the atoms are crystal- 
lized in the most dense isotropic structure, namely the 
FCC lattice. In order to stabilize other structures, sev- 
eral phenomenological potentials have been introduced, 
often without a microscopic justification. For example, 
the SW potential or the Tersoff potential stabilize the di- 
amond structure at low temperatures. These potentials 
have been used with success to calculate properties of Si 
clusters 13| and amorphous Si crystals 
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In this paper, we use the LJ potential to study the 
melting of rare-gas crystals which have the FCC lattice 
structure at low temperatures. It is amazing that such a 
simple question was not studied with precision so far in 
spite of an abundance of experimental data on rare gas 
such as Ar, Xe, Ne and Kr. Most of the melting studies 
concerning rare gas were done in particular cases: small 
clusters [14(, adlayers on a substrate, etc. The main rea- 
son to avoid to study the bulk melting may be due to 
some technical difficulties such as periodic boundary con- 
ditions, volume expansion with temperature, etc. Previ- 
ous Monte Carlo (MC) studies of bulk melting have been 
carried out with LJ potential but emphasize was put on 
the melting mechanism rather than on the precise melt- 
ing temperature in real materials 0, Q • 

The purpose of this paper is therefore to test whether 
or not the experimental T m can be reproduced by MC 
simulation using the values of the LJ parameters deduced 
for rare gas in the gaseous state long time ago[15|. We 
will show here that by appropriate choices of technical 
procedures, we are able to obtain melting temperature 
for various rare gas directly from our simulations, unlike 
previous simulations [lq - tl9 | which have had recourse to 
various means and some thermodynamic functions to de- 
duce it. We find in this work the melting temperatures 
for several rare gas higher than experimental values. A 
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revision of the values of LJ parameters widely used in 
the literature for more than 50 years should be made in 
order to better describe the solid state of rare gas. Note 
that in a recent work [20| , a hypothetical thermodynamic 
integration path is used to find the relative free energies 
of the solid and liquid phases, for various system sizes, 
at constant cutoff radius, in an attempt to explain the 
overestimate of the melting temperature with the LJ po- 
tential. However, due to various approximations, several 
results were not physically clear, in particular why the 
melting temperature oscillates with increasing cutoff dis- 
tance. 

Section |TI] is devoted to a description of the model and 
our MC technique. The results are shown and discussed 
in Sect. IIIII Concluding remarks are given in Sect. IIVI 
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Figure 1. (Color online) The Lennard- Jones potential for Ar 
with the parameters listed in Tabled] 



II. MODEL AND METHOD 



Table I. Lenard- Jones parameters 



The LJ potential is given by 
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with 



(1) 



(2) 



where e , a, b and a are constants, r»j denotes the dis- 
tance between two atoms at and Vj . Here a = 1 
and 6=1. We list for convenience the values of the 
above constants for various rare gas and their melting 
temperatures, using the data from Ref. Q . Note that the 
constants e and a are those deduced with some approx- 
imations using experimental data of the gaseous phase 
of rare gas which have been described in the pioneer pa- 
per of BernardesJUj]. The listed values of the constants 
should be therefore viewed as approximate values. To our 
knowledge, there were not papers using these constants 
to verify if the melting temperature is correctly obtained, 
either by theoretical calculations or by simulations. The 
absence of works dealing with this point has motivated 
our present work. 

In one of our previous works Q, we have used the 
LJ potential to investigate the mechanism which initi- 
ates the melting. We have counted the number of defects 
created with increasing temperature (T) and found that 
the melting occurs when these defects interact with each 
other forming a kind of chains of defects which break 
the solid periodic state at a temperature well below the 
melting point. In that work, we did not use the con- 
stants of the LJ potential for any specific material. We 
used instead the reduced temperature and dimensionless 
parameters. It was not our purpose to clarify the value 
of melting temperature for a given kind of crystal. It is 
now the time to verify those constants to see whether or 
not they yield the correct value of T m for the melting of 
the solid phase of rare gas. 



Element 


e (eV) 


<r(A) 


T m (K) 


Ar 


0.01042332126 


3.40 


84 


Kr 


0.01404339691 


3.65 


117 


Ne 


0.00312075487 


4.74 


24 


Xe 


0.01997283116 


3.98 


161 


He 


0.00087381136 


2.56 





Reference 



Let us describe our MC technique in the study of crys- 
tal melting. In our simulations we take a system of FCC 
lattice of size N — 4xLxLxL where L = 4, 5, 6, 8 
and 10 namely a system of N = 256 , 500, 1024, 2048 
and 4000 atoms. Unlike spin systems on lattice where 
the system volume does not change with T, a system of 
vibrating atoms expands its volume with increasing T. 
We use the method described in Ref. [2l[ but with a 
modification that the displacement of an atom i from its 
current position is taken in a sphere of very small ra- 
dius r with two arbitrary angles (9,<f>). To respect the 
uniform spatial distribution around r^, we use the Jaco- 
bian, namely we take an arbitrary cos(0) between -1 and 
+1, instead of an arbitrary 6. <f> is taken at random be- 
tween and 27T. Note that we first move all atoms to 
new positions, and change the system volume for a small 
arbitrary amount, then calculate the new system energy. 
We tune the magnitude of displacement (r) of the atoms 
and the volume variation in order to have an accepta- 
tion rate of displacement between 30% and 50% . This 
criterion is empiric but it is frequently used in MC sim- 
ulation. A MC step consists of moving all atoms, each 
with an arbitrary displacement, and changing the system 
volume once. The transition probability to the new state 
is exp(- 



fefr) where 



W — P (V new - V old ) + 2 (Unew - 
Nk B Tln 
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where P is the pressure which is set to zero here (constant 
pressure), V a id and V new are old and new system volumes, 
U id and Unew old and new system energies. 

One of the most difficult tasks in melting simulations 
is the application of the periodic boundary conditions. 
This is due to random positions of atoms at the crys- 
tal boundaries. We have taken the following actions to 
shorten, without loosing accuracy, our simulation CPU 
time: 

(i) each atom has a list of neighbors up to a distance 
Td longer than the potential cutoff distance r c . To es- 
tablish the list for the first time, we have to calculate all 
distances and we select neighbors at r < d. It takes time. 
The fact to choose rj, > r c is to ensure that for small dis- 
placements, some neighbors with r c < r < can enter 
inside the cutoff sphere without the need to reestablish 
the list of neighbors 

(ii) when an atom moves outside a sphere of radius r\ 
around its equilibrium position, its list of neighbors is 
updated by recalculating all pairs, neighbors up to fill 
the new list. r\ is chosen to be equal to 20% of ro the 
nearest-neighbor distance at equilibrium (T = 0) 

(iii) periodic boundary conditions are applied by trans- 
lating the system in all directions in a manner that avoids 
mismatches at the boundaries. 

The last step is the most important technically. We 
show in Fig. [2]a snapshot of the system and its translated 
atoms by periodic boundary conditions. We see that if 
they are of the same color, one cannot distinguish the 
system boundaries. 

In the following, we will show results with a cutoff dis- 
tance r c equal to twice £, the lattice constant of the FCC 
cell, namely r c — 2£ — 2roV2 where ro is the nearest- 
neighbor distance. As seen in section UlID I this value of 
r c is large enough to ensure a correct value of the melting 
temperature. 





Figure 2. (Color online) Upper: Snapshot of the system (ma- 
genta) with its translated atoms (gray) for periodic boundary 
conditions, Lower: The same snapshot of the system with its 
translated atoms but with the same color (online): bound- 
aries are undistinguishable. Argon crystal at T = 95 K with 
N = 500 atoms. 



to equilibrium depends on various MC parameters such 
as value of displacement magnitude r, volume variation <5 
etc. So, nothing can replace an observation of the time- 
dependence of physical parameters during the simulation, 
although this consumes a huge computer memory. 



III. MELTING OF RARE-GAS CRYSTALS 

A. The case of Ar 

Let us show first the results for Ar obtained by using 
the values listed in Table HI We will discuss next the mod- 
ification necessary for obtaining the result in agreement 
with experiment. In order to take a correct average of 
physical quantities, we record spontaneous values of all 
physical quantities during each MC run. We have to go to 
several millions of MC steps before observing statistical 
fluctuations around equilibrium. We show an example 
of the energy per atom E versus MC time in Fig. [3] for 
iV = 256 atoms at two temperatures T — 92 K and 94 
K. At T — 92 K, the system is still in the solid phase. 
Its energy is stabilized after about one million MC steps 
per atom. However, at T — 94 K, the system melts after 
three millions of MC steps per atom: E is stabilized in 
the liquid phase only after such a long MC time. It is 
very important to emphasize that the convergence time 
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Figure 3. (Color online) Energy per atom E versus number of 
MC steps per atom t for an Argon crystal at T = 92 K (red) 
and T = 94 K (green) with N = 256 atoms. 

We show in Fig. 0]the energy per atom E versus T in 
the case of Ar with N = 256 atoms, using the parameters 
given in Table HI We observe here that the melting occurs 
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at T m ~ 93 K with a large latent heat. This value of T m 
is rather far from the experimental data given in Table [I] 
We return to this point later. We show now the snapshots 
of the system for different temperatures in Fig. [5] below 
and above the melting. As seen, the system just starts 
to be spatially disordered at 92 K, and is well in the 
disordered phase (liquid) at 100 K. 
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Figure 4. (Color online) Energy per atom E versus tempera- 
ture T for an Argon crystal of 256 atoms, obtained with the 
parameters listed in Table [I] One observes that the melting 
temperature is T m ~ 93 K while the experimental value is 84 
K. Run of 50 millions MC steps/atom. See text for comments. 
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Figure 6. (Color online) Radial distribution of Ar for different 
temperatures: below melting T = 74 K (discontinued red) 
and 92 K (discontinued green), above melting T = 94 K (solid 
blue) and 100 K (solid magenta), with N = 256 atoms. 
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Figure 7. (Color online)Evolution of the lattice constant 
against temperature T. As we can see, there is a jump of 
this quantity which occurs at the same temperature of the 
jump of the energy in Fig. [4] Argon crystal, N = 256. 



Figure 5. (Color online)Instantaneous pictures of the super- 
cell of Ar for different temperatures, respectively from left to 
right : 75 K (far below T m ), 92 K (close to T m ) and 100 K 
(above T m ). 

The radial distribution for Ar is shown in Fig. [6] at 
several temperatures. One sees clearly the peaks at first, 
second and third neighbor distances for T < T m indi- 
cating the crystalline phase, while there is no such clear 
distinction for T > T m where a liquid phase is set in. 

We display in Fig. [7] the curve which shows the dilata- 
tion of the simulation box against the temperature. 

Let us discuss now the difference between our value 
of T m (N = 256) = 93 K with the experimental value 
T m (exp)=84 K. Clearly, the parameters given in Table 
U which have been deduced with experimental data for 
gaseous Ar do not describe well the melting of solid Ar. 
Note however that the above value of T m is for N = 256 
atoms. If we consider larger samples, T m should increase 
further. 

In order to modify correctly the LJ parameters while 
respecting the known properties of Ar, we use the listed 



values of LJ parameters to estimate the size effect, 
namely to find the value of T m (N = oo) for an infinite 
crystal. Once this task is done, we then we look for the 
parameters (e, a) which give correctly the experimental 
T m (exp). The effect of the cutoff distance will be exam- 
ined. 



B. Size effect 

The size effect is an important fact when we work with 
numerical simulations because of the finite size of the 
simulation cell. In second-order phase transitions, finite- 
size scaling laws allow us to calculate the critical expo- 
nents by varying the system size [22T[25l | . The correlation 
length is infinite at the critical point in the thermody- 
namic limit. However, in first-order transitions, the cor- 
relation length is finite at the transition point, the two 
phases coexist [2ol [27j and the energy is discontinuous. 
The correlation is finite at the transition temperature 
means that when the system size exceeds the correlation 
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length the transition temperature does not depend any- 
more on the system size. We need not to go to the infinite 
limit to find saturated transition temperature. Note that 
if the size of the system is too small, a first-order phase 
transition can appear as a second-order phase transition 
if the correlation length exceeds the system size. As we 
can see in FigjSJ fortunately when the size of the simula- 
tion cell is larger than 5 lattice constants, melting tem- 
perature T m is bounded by 102 K. It means that the 
correlation length at the melting temperature is smaller 
than 5 lattice spacings. 



Table II. Experimental (Exp.) and theoretical (Th.) values 
of NN distance ro , cohesive energy uo and bulk modulus 73cQ 
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Figure 8. (Color online)Energy per atom versus T with var- 
ious system sizes from TV = 256 (4 FCC lattice cells) to 
N = 4000 (10 FCC lattice cells). The arrows indicate T m 
for the smallest and largest sizes. See text for comments. 

The saturated melting temperature obtained by our 
MC simulations is thus 18 K (or 20%) higher than the ex- 
perimental one. This disagreement come from the values 
of the Bernardes parameters used for the LJ potential. 
The reason is that these values have not been calculated 
by fitting with the melting temperature but they were 
calculated with a low-density gas using the second virial 
coefficient [H]]. The Bernardes parameters are therefore 
questionable as we can see in Refs. 28-30]. Taken from 
those articles, the values of the nearest-neighbor distance, 
the cohesive energy and the bulk modulus calculated us- 
ing the Bernardes parameters are listed in Table [TT1 We 
can see the differences between experimental values and 
theoretical ones for all cases Ne, Ar, Kr and Xe. The 
Bernardes parameters also yield the high melting tem- 
perature found in our simulations with respect to the 
experimental one. 

At this stage, it is interesting to note that Molecular 
Dynamics (MD) simulation of melting of a perfect crystal 
with periodic boundary conditions produce superheating. 
There is an empirical rule which states that the melting 
temperature of a crystal without any surface and any de- 
fects, is 20% higher than the true thermodynamic melt- 
ing temperature T m . However in MC simulations, defects 
and dislocations are naturally created in the crystal by 
means of random numbers used in every MC step for 
atom displacements. Thus, the superheatin g sh ould not 
exist. Agrawal et al. have shown in Ref. 31( that for 
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Ar, with MC simulation, T m is about 15% higher than 
the experimental value. With our results for Ar, we find 
an increase of about 20%. For Kr, we find an increase 
of 13% for N = 256 atoms after 23 x 10 6 MC steps 
per atom. This increase is more important if we consider 
larger sizes as seen in the case of Ar. As said, the high 
values of T m in MC simulations are not due to the super- 
heating as in MD simulations. Rather, we believe that 
these high values are due to the inaccuracy of the listed 
Bernardes parameters. We will propose a modification 
in the following 



C. Modification of Lennard- Jones parameters 

In order to reduce the T m value, we propose now to 
modify the value of e, the prefactor of the LJ potential, 
and the coefficient a. 

Note that in papers dealing with melting in other ma- 
terials by means of MD calculations or MC simulations, 
there have been several propositions to modify constants 
appearing in potentials in order to obtain a correct agree- 
ment with experimental value of T m . This is because 
these constants are often deduced from experimental data 
which are not valid for the whole temperature range. 
Among these papers, we can mention the case of melting 
of Si crystal studied by MD [32], [HJ and MC simulations 
[3lL l34l ] , using the Tersoff potential [ll[ . Our proposition 
to modify some constants of the LJ potential when ap- 
plied to a rare-gas crystal is certainly a necessity in order 
to reproduce the experimental T m . 

We have done simulations with different pairs of (e, a). 
It turned out that T m depends essentially on e. There is 
however an optimal value of a which is a — 3.44 A corre- 
sponding to the experimental nearest-neighbor distance 
r = 3.75 A of Ar (cf. Tabled. We show in Fig. |9]the 
curves obtained for two selected values of e =0.008767853 
and 0.008951794 which give respectively T m = 83 K and 
86 K. These values oiT m are in agreement with the exper- 
imental value 84 K within statistical errors. Note that 
the modified e is about 15% smaller than the original 
Bernardes value e = 0.01042332126. 
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Figure 9. (Color online) Energy per atom versus T for two 
selected values of e, with a — 3.44 A, for an Ar crystal with 
N = 500. The arrows indicate T m for the two indicated values 
of e. See text for comments. 



D. Effect of cutoff distance 



At this stage, a natural question we ask ourselves is 
"what is the effect of the cutoff distance r c ?". We know 
that for a long-range interaction, the longer the interac- 
tion range is the lower the energy becomes. As a con- 
sequence, the melting transition is higher. However, as 
r c increases, the contribution of neglected neighbors be- 
comes smaller. From a certain value of r c , T m does not 
vary significantly. This is observed in Fig. [10] where T m 
is saturated for r c > 21, i. e. r c > 10.6 A. All the results 
shown above for r c are valid in the discussion of the size 
effect and the modification of e and a. 



E. The case of other rare gases 

In order to show that our algorithm works well with 
other rare gas, we have plotted the curve of energy versus 
temperature, obtained for Krypton and Xenon in Fig. 
[TT1 Again here, we see that T m , even for a small size, 
is already higher than the experimental value for each 
crystal. We think that the Bernardes parameters for Kr 
and Xe should be modified to get an agreement with 
experiments as what proposed above for Ar. 
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Figure 10. (Color online) Energy per atom versus T with 
various values of r c for Ar crystal with N = 500. The left 
arrow indicates T m for r c — 1.51 and the right arrow indicates 
T m for r c = 21 ~ 10.6 Aand 3£ ~ 15.9 A. Note that t is the 
FCC lattice constant which is equal to roy/2 where ro = 3.75 
A is the NN distance. See text for comments. 



Figure 11. (Color online) Upper: Energy versus temperature 
for a Krypton perfect crystal with N = 256 atoms. One 
observes that T m = 132 K while the experimental value is 
117 K. Lower: Energy versus temperature for a Xenon perfect 
crystal with N = 500 atoms. One observes that T m = 191 
K while the experimental value is 161 K. These curves have 
been obtained with the Bernardes values of parameters. 



IV. CONCLUSION 

We have shown in this paper results of the melting 
temperature for rare gas, by performing extensive MC 
simulations with the LJ potential. We have obtained di- 
rectly from our simulations physical quantities such as 
internal energy, lattice constant and radial distribution 
as functions of temperature. We have shown that melting 
occurs with a large latent heat and a jump in the lattice 
constant. Effects of system size and cutoff distance have 
been investigated. Let us emphasize that the theoretical 
LJ parameters widely used in the literature fl5j yield a 
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melting temperature higher than the experimental one as 
seen above, from ~ 15% for Kr to ~ 20% for Ar. This 
is not a surprise because those LJ parameters already 
yield theoretical NN distance, cohesive energy and bulk 
modulus different from corresponding experimental ones 
(see Table HIl). We have demonstrated that, in order to 
reduce the melting temperature to fit with experiments, 
it is necessary to modify the original Bernardes LJ pa- 
rameter e in such a way to reduce the energy at T = 0. 



The effect of er, within possible values of NN distance, 
is very small on T m . A good agreement on T m between 
experiments and simulations for Ar is obtained with the 
modified values given in Fig. O For other rare gas such 
as Kr and Xe, we have to proceed to a modification of 
their Barnardes parameters in a manner similar to that 
done above for Ar if we want to get a good agreement 
with experiments. 
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